Error thresholds for self- and cross-specific enzymatic 

replication 



Bcncdikt Obcrmaycr a , Erwin Frey a '* 

a Arnold- Sommerf eld- Center for Theoretical Physics and Center for NanoScience, 
Ludwig-Maximilians-Universitdt Miinchen, Theresienstr. 37, 80333 Miinchen, Germany 



Abstract 

The information content of a non-enzymatic self-replicator is limited by Eigen's 
error threshold. Presumably, enzymatic replication can maintain higher com- 
plexity, but in a competitive environment such a replicator is faced with two 
problems related to its twofold role as enzyme and substrate: as enzyme, it 
should replicate itself rather than wastefully copy non-functional substrates, and 
as substrate it should preferably be replicated by superior enzymes instead of 
less-efficient mutants. Because specific recognition can enforce these propen- 
sities, we thoroughly analyze an idealized quasispecies model for enzymatic 
replication, with replication rates that are either a decreasing (self-specific) or 
increasing (cross-specific) function of the Hamming distance between the recog- 
nition or "tag" sequences of enzyme and substrate. We find that very weak 
self-specificity suffices to localize a population about a master sequence and 
thus to preserve its information, while simultaneous localization about comple- 
mentary sequences in the cross-specific case is more challenging. A surprising 
result is that stronger specificity constraints allow longer recognition sequences, 
because the populations are better localized. Extrapolating from experimen- 
tal data, we obtain rough quantitative estimates for the maximal length of the 
recognition or tag sequence that can be used to reliably discriminate appropriate 
and infcasiblc enzymes and substrates, respectively. 

Keywords: origin of life, quasispecies theory, enzymatic self-replication, 
higher-order catalysis 



1. Introduction 



The acclaimed experimental finding (jCechl . llQQOh that RNA not only stores 
genetic information but also provides catalytic function has inspired the RNA 
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world scenario ( Gilbert , 19861) . a hypothesis for the starting point of Dar- 
winian evolution at the origin of life through self-replication of RNA-like polynu- 
cleotides. Substantial progress over the last decades demonstrates t he capability 

of RNA enzymes ( r ibozymes) to catalyze diverse chemical reactions (jDoudna and Cecb 
23 



2002t iLillevl . l2005t iJovcel 120071 ) , amon g them the polymerization of as many as 



20 nu cleotides to a template molecule ((Johnston et all 120011 ; iZaher and Unraul . 



20071) , and even replication through templa te-directed ligation react ions involv- 
ing short RNA oligomers as building blocks ( Lincoln and Joyce . l2009h . However, 
the currently known systems are not yet capable of Darwinian evolution, lacking 
either the ability to replicate molecules as long and complex as themselves or 
to introduce heritable variation. 

For theorists, the focus has mainly been on whether the information content 
of a self-replicating molecule can be maintained in the presence of replication 



error s, usually employing Eigcn's well-known quasispccics theory (jEigen et al 



19891) for the self-replication of L-nucleotide sequences. In this model, replica- 
tion errors occur with an error probability \i per single nucleotide, and the repli- 
cation rates are taken as functions only of the template sequence. This mapping 
from genotype to replication rate (or fitness) presents a considerable challenge 
for evolutionary models, and even generic structural features of such fitness 
landscapes, such as the extent of ruggedness and epistasis, are under debate. 



Valua ble insight comes from experiments on mutagenized ribozymes (|Kun et al 



20051) a nd from computer simulations using RNA secondary structure as phe 
notype ( Huvnen et al. . 19961 : Takeuchi et al. . 2005 ). which indicate a significant 
degree of neutrality around an optimal sequence, reducing the deleterious ef- 
fects of mutations. Other theoretical studies used a large number of different 
idealize d fitness landscapes, e.g., with a single peak at one fittest "master" se- 
quence (ISwetina and Schuster . 19821 Leuthausser , 1986 ; Schuster and Swetina . 



Woodcock and HiggsHl996HGallucc 



±1 
al., 



Hermisson et all 120021; IPeliti 



1988; 



2002: 



1997 



20091) , with some rather generic re- 



. .Saakian and Hul. 12006 ; Saaki an et all . 
suits (jWiehel . Il997t IJain and Kruei l2005l) ~ the population in sequence space is 
characterized by a broad mutant distribution (a quasispecies) localized about 
the master sequence for mutation probabilities smaller than a critical value /x c 
(the error threshold), while it consists of random sequences (it is delocalized) 
for larger values. Because the error threshold [i c ~ l/L is usually inversely 
proportional to sequence length, the problem arises whether the maximally sus- 
tainable comple xity of a self-replicator su ffices to perform the complex task of 



19781) . 



self-replication (jEigen and Schuster . 

In a prebiotic context, it is important to emphasize that using a fitness 
landscape where rates depend only on the template sequence pertains to non- 
enzymatic rather than enzymatic replication, because in the latter case the 
replication rates also depend on the concentrations and the characteristics of 
involved enzymes. For RNA, the potential for non-enzymatic replication is 
questionable, given that template-directed polymerizatio n or ligation seems lim 



ited to short molecules with rather specific sequences (|von Kiedrowskil 11986 



Acevedo and Orgell . Il987l: IWu and Orgell . 119921 lOrgell . 12004) . Moreover, experi- 
mental and simulation studies demonstrate a strongly disadvantageous tendency 
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for elongating side-reactions at the cost of repli cation (IFernando et al 



Hence, although there have been speculations ( Pace and Marshl Il985 ) , it re- 



2007). 



mains unclear how a single more complex RNA should literally copy itself ( Jovcd . 



2007t ISzostak et all 12001 ) 



Enzymatic replication is more plausible (|OrgelLll992ft . but raises the question 
whether high replication efficiency (high fitness) is a property of the substrate 
or the enzyme. In the latter case, a superior rcplicase does not enjoy a selective 
advantage, because it replicates non-functional mutant templates just as well as 
itself, while it is not guaranteed in the former case that a superior template is 
functional at all. Likewise, mutations generate substrate s that are replicated les s 
efficiently, but they also produce less-efficient enzymes (jMavnard Smith! . Il979h . 
thus affecting the replication rates of all potential substrates. On theoretical 
grounds, one should expect that a superior replicator is a good enzyme and a 
good substrate at the same time. As enzyme, it should therefore replicate only 
functional substrates, and as substrate, it should be replicated preferably by 
efficient enzymes. It has long been realized that in a competitive environment 
these propensities are crucial for the emergence, improvement and perpetua- 
tion of replicase a ctivity, which i s essentially an altruistic trait that is not by 



itself selected for ( Michod . 19831 ). The commonly proposed solution to ren- 



der enzymatic replication evolutionarily stable is to impose a form of group 
selection, e.g., via compart mentalization in vesicles, in order to keep s imilar 



molecules closely toge ther (jSzathmarv and Demeterl . Il987t lAlves et al.l . 12001 



Fontanari et all 120061 ) . even though this requ ires a simultaneous and coordi- 
nated emergence of replicators and protocells (jSzostak et all l200ll ). However, 



known ribozymes a ct with moderately or even strongly substrate-specific effi- 
ciency ([Jovcd . 120071 ) , such that specific recognition may also have a significant 
influence. Since unspecific reactions require sophisticated substrate-binding 
proper ties that could well ha ve been a rather late invention in prebiotic evo- 
lution ( Johnston et al. ■ l200lh . it seems natural to assume that replication effi- 
ciency should depend strongly on the interaction between enzyme and substrate. 

In this paper, we analyze a general model of enzymatic replication account- 
ing both for varying degrees of specificity and the characteristically broad quasi- 
species distributions in order to address the consequences for the error threshold . 
Simil ar to models for the evolution of regulatory DNA motifs ([Gerland and Hwa , 
2002), we assume that specificity depends on the quality of binding to some 



recognition or tag sites (jWeiner and Maizeld . [l987T) . Idealizing this condition, 
we use replication rates that depend on the Hamming distance between these 
sequence regions of enzyme and substrate via a decreasing (self-specific) or in- 
creasing (cross-specific) function. After formulating the model and discussing 
our analytical approach in a methods section, we show for these two scenarios 
results from stochastic simulations and numerical solutions of deterministic rate 
equations. Using our analytical toolbox, we discuss the resulting localization 
conditions, error thresholds and the phase diagram. In our conclusions, we use 
experimental values for polymerization rates to obtain simple estimates for the 
maximum number of nucleotides that can be used for recognition. 
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2. Model 



In the framework of quasispecies theory, each molecule is characterized by its 



sequence Si = (a^ . . . cr^') of L binary nucleotides a\ L> € {0, 1}. In an infinitely 
large popu lation, its concentr ation Xi evolves according to the deterministic rate 



(0 



equations (jEig en et all Il989h 



(1) 



Here, My 



is the mutation probability between sequences 



Si and Sj with Hamming distance = yj £ jcr^ — er^l, where fi is the er- 
ror probability per single nucleotide (usually called "mutation rate"), and the 
replication rate Ri of sequence Si is given by: 



JJ) i 



R, 



(2) 



Whereas the non-enzymatic rate A4 depends only on the genotype Si, the second 
term implies frequency-dependent selection and makes our model intrinsically 
nonlinear. It encodes the catalytic interactions of two molecules: By measures 
how well Sj catalyzes the replication of Si. The second term in Eq. (JTJ en- 
sures the normalization y\ Xj = 1, and a degradation term —DiXi therefore 
drops out of Eq. (JTJ since we assume that the decay rate Di = D is sequence- 
independent for simplicity. Note that actual RNA sequences replicate via a 
complementary intermediate, while our idealized model assumes direct replica- 
tion in a single step. It has been shown that these two ap proaches are es sentially 
equivalent for the symmetric situations considered here (jStadlerl . 119911) . 

To capture the pertinent features of a situation where replicase enzymes pre- 
fer to replicate themselves instead of their competitors, we assume that the qual- 
ity of specific recognition influences catalytic rates more strongly than the actual 
genotypes of enzyme and substrate. Hence, we effectively only model the recog- 
nition reg ions of riboz ymes, which are often clearly separated from the catalytic 
domains (|Lillevi . 120051) . Of course, the proper function of the latter region is in- 
dispensable, but in our idealized model we neglect the influence of mutations: as 
we have shown previously in a simple model, their effect on the error threshold 
is largely independent from the more interes t ing co nsequences of mutations in 
the recognition region ( Obermaver and Frev . 20091 ). The sequence length L is 
thus restricted to the number of nucleotides that take part in recognition . Prob - 
ably mediated via specific base-pairing interactions (jDoudna and Cechl . 120021 ). 
the quality of recognition can be taken as function of the number of mismatches 
between the binding sites of enzyme and substrate, and we let the catalytic ma- 
trix Bij therefore depend via a specificity function f(d) only on the Hamming 
distance dij between enzyme and substrate. Furth er, because rate enhance - 



ments through ribozyme catalysis can be substantial (jDoudna and Cechl . 120021 ) . 
such that non-enzymatic replication rates are comparably small (if nonzero at 



4 



all), we neglect their genotype dependence altogether and choose a flat fitness 
landscape for Af. 

A = a, By = (3) 

Because we are interested in the stationary state, the parameter a (if nonzero) 
merely sets the time scale while (3 measures the selection strength. 

We will first analyze a scenario for self-specific replication, where replication 
rates increase with similarity of enzyme and substrate. To have the degree of 
specificity explicitly tunable via a parameter p, we use the specificity function 

/.(d) = (1 - d/Lf. (4) 

As a contrasting example, we will then analyze a similarly defined specificity 
function where replication rates increase with complementarity between enzyme 
and substrate: 

f c (d) = (d/Lf. (5) 

The numerical and analytical approach presented in the next section can readily 
be applied to other functional forms of the specificity function. 



3. Methods 



For reasonably large sequence length, the full 2 -dimensional system Eq. ([TJ 
can only be analyzed using stochastic simulations in a finite population of N 
sequenc es. Here, we e mploy the straightforward stochastic simulation algorithm 

fk copies, 



sequenc es. Here, we e mploy the straigntiorward stochastic simulation alj 
used by IWilke et al.l ( 200 ll ) . At time t each sequence Sk , present in Nk 



has a probability po,fc = ^k/J2i A^(l + Ri) to be copied without mutations into 
the population at time t + 1, and a probability p mu t,jk — MjkRkNk/ Y^i + 
Ri) to be selected and mutated into sequence Sj. Following initialization, our 
observables of interest are measured by averaging over time after reaching a 
stationary state. 

In order to derive analytical results, we exploit that the stationary states of 
Eq. ([I} are localized about a particular "master" sequence S 1 * , which is necessary 
to preserve its information content. In contrast, derealization indicates that 
such a sequence S* cannot be maintained due to replication errors. Since the 
possibility of non-trivial dynamics such as periodic orbits cannot be excluded 



for general replicator- mutator equations like Eq. (p} (jStadler et all 1 1995T) , there 
might also be other reasons for the absence of localization. However, we did 
not find any signs of complex dynamical behavior in our simulations, and one 
can easily convince oneself that the delocalizcd state, where all sequences have 
the same concentration and therefore identical replication rates, can lead to lo- 
calization: since replication rates are essentially proportional to concentration, 
stochastic concentration fluctuations imply higher rates. Unless hindered by 
excessive mutations, this can induce a transition to a state localized about some 
randoml y chosen master sequence, w hich then also has the highest replication 



rate (see Obermaver and Frev ( 20091 ) for a visualization). Such "fixation" events 



are very similar to the phenomenon of consensus formation, e.g., in language 
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dynamics (<Blvthd . [2009h . Note that our idealized replication rates do not pre- 



determine any specific master sequence for localization. This symmetry would 
be broken in a full model where replication rates depend on the full genotypes 
of enzyme and substrate (and not just the Hamming distance between their 
recognition regions). 

Given localization about S 1 *, we can significantly reduce the dimensionality 
of Eq. ([T]) by lumping all sequences Si with a Hamming distance k to S* together 
into "error class" k. Witho ut loss of generality, we assume that S* = (00 . . . 0) . 
This well-known procedure ( Schuster and Swetina . 1988 ; Woodcock and Higgsl . 



19961 ) permits to formulate reduced rate equations formally equivalent to Eq. |T]) 
in terms of new variables Xk denoting the concentration of error class k in the 
population: 

x k = y^^mkjjajSji + bjiXi\xj - x k ^7" A" + b ji x l ]x j . (6) 

ji ji 

Here and in the following, we use lowercase letters for all variables in the reduced 
system. For our model Eq. (|3]), the non-enzymatic rates are given by eij = a. 
The accordingly reduced mutation matrix and the catalytic matrix depend only 
on the Hamming distance between pairs of sequences in different error classes 
(measured with respect to the master sequence), which allows us to combinato- 
rially assess all possibilities for their relative distance. The total probability of 
distributing 0<k — j + 2£<L mutations to move a sequence from error class 
j into error class k has been derived previously as ( Woodcock and Higgsl 1 1996T) 



The replication rate Eq. ([2]) also depends on the frequency of each sequence in 
each error class. With the homogeneity assumption that all (^) sequences in 
class j are equally populated, this complication can be resolved, and the reduced 
matrix bij reads analogously 

Numerical solutions to the (L + l)-dimensional rate equations given in Eq. ^ 
can now easily be found by means of standard algorithms. 

A more detailed understanding of these solutions can be obtained from the 
population distribution's normalized first moment a = (k) / L = (a), which as 
the mean Hamming distance (k) = kxk to the master S* = (00 ... 0) char- 
acterizes the width of the distribution and measures the mean value (a) of each 
sequence's binary nucleotides. Writing down an equation for the first moment 
of Eq. (|()]) requires a hierarchy of expressions for higher moments, which can 
be truncated by means of a moment closure technique. For the self-specific 
case, we assume that the stationary Hamming distance distribution is approx- 
imately binomial, Xk ~ (^)a fc (l — a) L_fc , because this reproduces the expected 
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distribution in the limits a — > (complete localization about one sequence) and 
a — > 1/2 (the delocalized state, where the binary nucleotides are random num- 
bers). Moreover, it solves th e rate equations Eg. Q exa ctly for linear fitness 
landscapes without epistasis ( Woodcock and Higgsl 1996) and for an ex tension 
of the quasispecics model to a game theory setting (jLassig et all 120031 ). With 



this binomial ansatz, a is the population distribution's only parameter, and for 
our model of the replication rates Eq. (j3)), it obeys the equation 

(1 - 2a){nL [a + £5(1 - 2a(l - a))} 

- a(l - o)(l - 2n)PS'{l - 2a(l - a))} = 0. (9) 



This equation for a, which is one of our main analytical results (see |Appcndix A.l| 
for a derivation), holds for any specificity function f(d), which enters via the 
function 



(10) 



An intuitive interpretation of this auxiliary function derives from recognizing 
that the quantity a + (3S(l — 2a(l — a)) measures the mean replication rate (or 
mean fitness) J^j RjXj of the population if a solves Eq. ([9]). Such solutions a(/x) 
can be obtained explicitly only in special cases where S(x) attains a simple form 
(see |Appendix A."2j ) , but we can easily solve for /j,(a) and invert graphically to 
obtain a bifurcation diagram: 



H(a) 



a + /3S(l - 2a(l - a)) 
/3a(l-a)5'(l-2a(l-a)) 



(11) 



For cross-specific replication with the specificity function Eq. ([5]) , we expect 
two equivalent subpopulations localized about complementary sequences, which 
corresponds to a superposition of binomial distributions: Xk ~ ^(^)[ ak {^ ~ 
a) L ~ k + (1 — a) fc a i_fc ] . To obtain an equation similar to Eq. ^ for their mean 
widths a, we cannot use the first moment of the reduced rate equations (it 
vanishes by construction, because the distribution is symmetric about a = 1/2), 
but use the second mo ment (Afc 2 ) = ^ k {k — (k)) 2 xt, leading to a lengthy 
result explicitly given in Appendix B.l| and a corresponding expression for the 
bifurcation diagram ( see Eos. (|B.5|) and (|B~7| ). 



4. Results and Discussion 

4-.1. Self-specific replication 

Our first scenario is concerned with self-specific replication, where replication 
rates increase with similarity of enzyme and substrate through the specificity 
function Eq. As argued in the preceding section, stochastic fluctuations 
increasing the frequency of one particular sequence also increase its replication 
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rate, such that the population can localize about this master sequence S*. Fig- 
ure [T] shows the stationary Hamming distance distribution Xk , which measures 
the concentration of sequences with k mutations relative to S 1 *, for different 
degrees of specificity from the linear case p = 1 to complete self-specificity 
p — > oo, where enzyme and substrate have to be identical. We compare results 
from stochastic simulations to numerical solutions of the reduced rate equa- 
tions, Eq. (|6|). The excellent agreement between simulation and deterministic 
theory justifies the homogeneity assumption made in symmetrizing the speci- 
ficity matrix (see Eq. ([5])). In the simulations shown in Fig. [TJ we initialized all 
sequences at a predetermined master sequence in order to avoid noise from the 
intrinsically stochastic "fixation" events, but we also tested other initial condi- 
tions with modest inhomogeneities, which were quickly "washed out" and did 
not give rise to measurable differences in the stationary state. 

The limit p — >• oo with f s (d) — > Sd,o, depicted in Fig.QTc), leads to a general- 
ize d Schlogl model of a uto-catalytic replication, which has been partly analyzed 
bv lStadler et al.J ( 1995). In this limit we can employ the well-known "error-tail" 



approximation JSchuster and Swetina . Il988l) : we define xq as the concentration 
of the master sequence, a + fix® its replication rate and (1 — fi) L the proba- 
bility not to have a mutation. All other sequences are lumped together in the 
error tail with concentration 1 — xq and replication rate a (the concentration of 
suitable replicase enzymes is so small that the frequency-dependent term in the 
replication rate does not contribute) . Neglecting back mutations from the error 
tail into xq (corresponding to the large-genome limit), we obtain the simple 
equation 

x Q = (a + f3x )(l - n) L xo - x f, (12) 

with f = (a + [3xq)xo + a(l — xo) the mean replication rate. In the sta- 
tionary state, we easily find that the delocalized state xo = is stable for 
all fi, while a branch of solutions with nonzero xq emerges f or (1 — n) L > 



2( ^/ a ( a + 0) — a) / through a disc o ntinuous trans i tion (s ee also lCampos et al 



(|2000h : lObermaver and Frevi (|2009t ) : IWagner et al.l (|201fj| ) for similar results in 



related models) . Whereas for the somewhat related sharply-peaked fitness land- 
scape, where only the master sequence has a hig her replication rate, the error 
threshold arises through a continuous bifurcation (jBaake and Wiehd . 11997ft . the 



discontinuity observed here expresses the qualitati vely different behavior we pre - 
viously termed "escalation of error catastrophe" ( Obermaver and Frev . 2009h : 



as the mutation rate grows, the proportion of fittest sequences, i.e., of the neces- 
sary replicase enzymes, is diminished and therefore their replication rate. This 
in turn reduces their concentration, until at the error threshold the concentra- 
tion of enzymes xq is not large enough to have them replicate with an efficiency 
sufficient for localization. 

Although Eq. (fT2"j) approximates the exact result for xo very well (see the 
dashed line in Fig. QIc)), it is valid only for L — > oo (because back mutations 
are neglected) and p — >• oo (because the replication rate of the error-tail is taken 
as concentration- independent). In order to gain a more general perspective, we 
use the analytical solutions of Eq. ([9]) for the mean Hamming distance to the 




Figure 1: Solutions for self-specific replication with f s {d) = (1 — d/L) p . Stationary Hamming 
distance distribution x^ as function of mutation rate fi from numerical solutions to the reduced 
rate equations Eq. $6$ (straight lines) and a stochastic simulation of the full system Eq. JTJ 
in a population of N = 10 4 sequences of length L = 32 (dots) for a = 1, f3 = 5 and (a) p = 1, 
(b) p = 5, and (c) p = oo (here, the dashed line shows the error-tail approximation for xo). 
The insets depict the average Hamming distance a = (k) /L to the master sequence. 
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master (the population distribution's first moment a = (k) /L) obtained from 
the bifurcation diagram Eq. A comparison between the exact solution a(a) 
resulting from the reduced rate equations Eq. ^ via numerical continuation^ 
and Eq. (fTTj) is shown in Fig. [5] for L = 8, a = j3 = 1 and different values of 
p in the specificity function Eq. (QJ. Recalling the symmetry of the original 
model, Eq. ([3]), namely that the population can localize about any sequence, 
the remaining reflection symmetry a — > 1 — a about the dclocalizcd solution 
a = 1/2 indicates that after symmetrization localization is only possible about 
the master sequence or its complement. The stable branches associated with 
localized solutions start for // = at a = (or a = 1), i.e., full localization 
about one sequence (or its complement), and higher mutation rates give rise to 
broader distributions with larger mean, until these localized regimes disappear 
at critical mutation rates jj, c denoted by circles. The delocalized solutions, on 
the other hand, gain stability at finite mutation rates jx c (denoted by crosses), 
and we find bistability for jl c < j.i < p c . A si milar situation is encountered in 
form ally related models of grammar evolution ( Nowak et all l200lt iKomarova , 
|2004 - 

Explicit expressions for the error threshold /i c = max /i(a) are available from 
Eq. (jTTJ) when the auxiliary function S{x) (defined in Eq. (|10[) ) has a simple form 
(see | Appendix A."2] ). For instance, /x c = /3/ (AaL + 2/3 (L + 1)) for p = 1, where 
specificity decreases linearly with the distance between enzyme and substrate. 
Altho ugh an error threshold is abs ent for linear fitness landscapes without epis- 
tasis (Wo odcock and Higgsl. Il996h . here the intrinsically nonlinear model gives 



a sharp transition even in this apparently similar case. For very strong speci- 
ficity p — > oo we find p c = W[P/(ea)]/(2L) using Lambert's W-function, and 
we recover the result p c = ln(/3/a)/(2L) that can also be obtained from the 
error-tail approximation in the limit f3 ^s> a. Further, the delocalized state 
a = 1/2 is the only solution of Eq. (J9j) in the complete absence of specificity 
(f s (d) = 1), supporting the intuition that unspecific replication docs not suffice 
to preferentially maintain the information of one particular sequence. Interest- 
ingly, taking p — > in Eq. (01 gives the finite even though exponentially small 
value p, c = (3/ (2 L+1 (a + /?)). This result implies that limited localization is pos- 
sible even for very weak specificity (if p = in the specificity function Eq. ^j, 
enzymes replicate everything except their exact complement, because always 
fs(L) = 0). 

From the bifurcation diagram obtained via Eq. (|11[) . we easily read off ex- 
act results for the value /2 C = m(1/2) where the delocalized state gains stabil- 
ity. E.g., for the genera lized Schlogl model p — > oo, we get p c = (3/(a2 L+1 + 



4/3) (jStadler et alj . 119951 ). This exponentially small yet finite value is consistent 
with our previous conclusion that the delocalized regime is stable for all values 
of p within the error-tail approximation, Eq. (|12p . which holds for L — > oo. Fur- 
ther, we find that the two critical values p c and p c are identical for p = 0, 1, 2. 
Recognizing that Fig. [5] describes a pitchfork bifurcation at a = 1/2 and fj. = fi c , 



AUTO software package available via http://indy.cs.concordia.ca/auto/. 
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we infer that the two critical mutation rates are equal (p c = p. c ) whenever the 
pitchfork is supercritical, whereas bistability between localized and delocalized 
states for intermediate mutation rates fl c < /i < fx c is possible in the subcritical 
case, leading to the discontinuous transition observed in Fig. [TJ The bistabil- 
ity regime vanishes as the curvature ^"(1/2) in the bifurcation diagram changes 
sign, which gives from Eq. an approximate expression for the corresponding 
critical value of 8: 



8* 



S ' 2 ^ 5(1/2, 



5"(l/2)-25'(l/2) 



(13) 



This general result can readily be evaluated for any specificity function entering 
the auxiliary function S(x) defined in Eq. (fpQ|) . It predicts bistability for all 
values 8 < 8* = a2 L ~ 1 (L — 2) if p — > oo in our choice Eq. (j4|) of the specificity 
function, and no bistability for weak specificity because the critical value of the 
coupling constant is negative (8* < 0) for small p < p m j n = 2 + Since 
the mean Hamming distance (k) = aL is discontinuous at the error threshold 
only in the bistability regime 8 < 8*{p), we argue that the observation of 
a discontinuous mean replication rate (mean fitness) a + 8S(l — 2a(l — a)) 
de pends on the specific ity of second order catalysis, which generalizes the result 
of IWaener etUl (l2010h . 



The results of the binomial closure approximation Eq. © are summarized 
in the phase diagram Fig. |3l where the two critical mutation rates [A c and fl c are 
shown as functions of the selection strength 8 and the specificity degree p for 
a = 1. The thick line denoted (3* indicates the boundary of the bistability regime 
/i c > /2 C . Fig. [5] demonstrates that the binomial approximation is quantitatively 
excellent in the supercritical situation 8 > 8* (in particular, it gives exact 
results for fi c and p, c ), and qualitatively correct otherwise, where the values 
for ^ c are somewhat underestimated: near the error threshold, the variance 
of the population distribution is considerably larger than that of a binomial. 
Nevertheless, we emphasize that the result ^ c = l n (/?/ Q; )/(2i) in the limiting 
case p — > oo agrees with the value obtained from the error-tail approximation 
for large 8. The remarkable performance of the binomial closure approximation 
can be appreciated in more detail from the projected phase diagrams shown in 
Fig. 12 

A noticeable feature of these phase diagrams is that the error threshold /x c 
increases for stronger specificity p (see Fig. HKb)), which implies that higher 
mutation rates can be tolerated, i.e., that longer sequences can be maintained. 
This seems at first counter-intuitive, because weaker specificity constraints on 
the recognition sequence should allow more mutational "freedom" . However, as 
shown in Fig. [5] and Fig. (TJ smaller values for p lead to much broader dis- 
tributions: mutants are still reasonably well replicated by master enzymes, 
but the master is only moderately well (but not quite as efficiently) repli- 
cated by the mutants. The resulting broadening of the distribution effectively 
redu ces the replication rate of the master and escalates the error catastro- 
phe ( Obermaver and Frev . 20091 ). Thus, the necessity for an enzymatic replica- 
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Figure 2: Bifurcation diagram of Eq. (O for the normalized mean Hamming distance a = 
(k) /L with L = 8, a = 1, /3 = 1 and different values of the specificity degree p. Thick lines 
indicate stable branches, dashed lines unstable branches and thin dotted lines the results of 
the binomial closure approximation Eq. J9j, which is barely visible for p = 0.1 and p = 1. 
Circles indicate critical mutation rates fi c where the localized regime vanishes, crosses show 
values jlc where the delocalized state changes stability (horizontal branches corresponding to 
delocalized states arc drawn slightly shifted for visualization). 




Figure 3: Phase diagram of localization regimes in the parameter space of mutation rate /i, 
selection strength f3 and specificity degree p (log-log-log scale) obtained from Eq. (|9} for L = 8 
and a = 1: below the upper plane fi c , a localized solution exists, while above the lower plane 
fi c the delocalized state becomes stable. Bistability (fi c > p, c ) is possible only for ft < /3*(p). 
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tor to discriminate not only between functional and non-functional substrates, 
but also between efficient and unproductive enzymes, is again emphasized. 

We finally want to remark on the case a = 0, i.e., no background level 
for the non-enzymatic replication rate. Most of the above results obtained 
from the binomial closure approximation can be simply evaluated for a = 
(note that then (3 sets the timescale and drops out), but the strong specificity 
limit p — > oo deserves extra attention. The error-tail approximation indicates 
that the error threshold vanishes {pL c — > 1), but from Eq. © we find the ex- 
act result ii c = fi c = 1/6 for L S> 1, i.e., a macroscopic yet finite value (see 
|Appcndix A.2[ ). This remarkab le result can be explained by rec alling that the 



traditional result p c w In r/L ( Eigen et al. . 19891 Wiehe . 1997 ) for the error 



threshold depends on the replication advantage r of the master relative to a 
possibly small but finite value for the mutants. In our case, rates are directly 
proportional to concentration, and because the master sequence has a concen- 
tration of order 1, while in an infinitely large population distant mutants have 
concentrations of order 2~ L , this relative advantage itself is of order 2 L , and 
cancels the length dependence of the error threshold. In the corresponding 
non-enzymatic case, results for so-called "truncation" fitness landscapes have 
lead to some debate about the ap plicability of the error threshold c o ncept 



in the presence of lethal mutations d Wilkc 



Takeuchi and Hogewed . 2007; Saakia n et al 



2005 
20091) 



Summers and Litwin. 2006 



Accordingly, we should cau- 
tiously note that our results for a = will probably be affected when accounting 
for the effects of finite populations and the full dependence of the replication 
rates on the genotypes of enzyme and substrate. 

4-2. Cross-specific replication 

To increase the information content of replicating systems beyond the limited 
complexit y of a single replicator, auto-catalytic reaction networks such as hy- 



percyclcs (lEigen and Schuster! . Il978t IStadler et all 119951) have been proposed, 



where different molecular species catalyze each other's replication in a possibly 
complex interaction graph. Only very little is known for these systems regard- 
ing the issue of reaction specificity and the cross-interactions of each species' 
mutant clouds. The simplest conceivable networks are 2-member cross-catalytic 
hypercycles. While the two members (which actually replicate via complemen- 
tary intermediates) need not be strictly complementary, we assume that they 
are sufficiently distinct that an idealized model requiring complementary recog- 
nition regions captures the essentials. This suggests to analyze the specificity 
function Eq. ([5]) where replication rates increase with Hamming distance be- 
tween enzyme and substrate. In this case, we expect the formation of two 
sub-populations localized about complementary sequences, each catalyzing the 
replication of the other. The main question to be answered is how specificity 
affects coexistence. In the standard Eigen model, quasispecies coexistence is 
prevented by competitive exclusion except i n degenerate cases, because the 
"fittest" individuals take over the population (jSwetina and SchusterLll982l) . In- 
teracti ons between sub-populations, e.g ., based on complementarity of binary 
traits ( de Oliveira and Fontanarl 20021 ). are known to enable coexistence. In 
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our case, we expect novel conditions for coexistence: even though each subpopu- 
lation depends on the presence of the other for efficient replication, it is unclear 
how the possibly broad mutant distributions influence each other. 

From numerical solutions to the reduced rate equations and simulation re- 
sults, where wc initialized the population split between the master sequence and 
its complement (see Fig. [5]), we find that localization is only possible for p > 1. 
This important result is confirmed through the binomial moment closure approx- 
imation (see Appendix B.2|, which has for p = 1 only the solution a = 1/2, 



corresponding to the delocalized state. Because the complementary distribu- 
tions overlap too much if catalytic rates increase only linearly with Hamming 
distance, the populations are not localized strongly enough to ensure coexis- 
tence. 

The bifurcation diagram /x(a) obtained from our approximation scheme (see 
Eq. (|B.7[0 is shown in comparison with the exact result from the reduced rate 
equations in Fig. [5] for different values of p. We obtain information about the 
parameter a from the variance (Afc 2 ) = ^2 k {k — (k)) 2 Xk, in the case of two 
complementary binomials given by (Afc 2 ) = L 2 /4 - L(L - l)a(l - a). Using 
this expression to infer a from the measured variance gives very good agree- 
ment between binomial closure approximation and exact numerical results es- 
pecially for small p, because now the first two moments are correct. The critical 
mutation rates fi c and fl c can be found from the bifurcation diagram, e.g., 
/i c = W[/3/(2ae)]/(2L) for p — > oo, which is identical to the corresponding re- 
sult for self-specific replication if we replace (3 — > (3/2. Hence, in this limit we 
obtain two clearly separated binomial distributions representing two equivalent 
and catalytically coupled populations: indeed, in the localized state the sum 
x k + x L-k from Fig. EJb) is equal to in the self-specific situation shown in 
Fig-HIc) once we replace f3 —> f3/2. The "coupling constant" is only half as large 
because only one half of the population is available as enzymes for the other. 
In particular, the negligible interaction between the respective mutant clouds 
allows one to emplo y the error-tail approximation assuming independe nt species 
and error tails as in ( Campos et al. . 2000t Obermaver and Frev . 2009h . Further, 



we find that the pitchfork bifurcation is always subcritical, i.e., that [i c > £l c . 
We summarize our main results by plotting projected phase diagrams of /x c and 
fl c as functions of p— 1 and /3 in Fig. [7] This confirms that both critical mutation 
rates /x c and p, c vanish linearly with p— 1, because n{a) en p — 1 as p — > 1, which 
gives the sharp bound p > 1 for coexistence of two populations. Finally, the 
case a — of zero non-enzymatic replication rate is similar to the self-specific 
case: the slight chance of distant mutants to find an appropriate enzyme gives 
an enormous replication advantage to the mainly populated master sequences 
and therefore macroscopic values for the error threshold (see |Appcndix B.3 for 
details). 

5. Conclusion 

Because enzymatic replication provokes the necessity for enzymes to favor 
functional substrates and for substrates to prefer efficient enzymes, we analyzed 
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Figure 4: Projected phase diagrams of localization regimes for self-specific replication with 
L = 8 and a = 1 obtained from Eq. J6j: (a) as function of with p = 10 fixed; (b) as function 
of p with f) = 10 fixed. Below fj, c (thick line), a localized solution exists, while above fi c (thick 
dashed line) the delocalizcd state becomes stable. The dotted lines denote the result for /z c 
obtained via the binomial closure approximation. Bistability (/i c > jl c ) is possible only for 
P<P*(p) (or p>p* 09)). 




Figure 5: Hamming distance distribution xj. as function of mutation rate /i as in Fig. [T] but 
for cross-specific replication with a = 1, f) = 10, L = 32, and (a) p = 5 or (b) p = oo. The 
insets show the variance (Afc 2 ) = *%2u[k — {k)) 2 Xk (in a symmetric population, (k) = L/2). 
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Figure 6: Bifurcation diagram of Eq. l(6j as in Fig. El but for cross-specific replication in 
terms of the population parameter a, which denotes the width of the two subpopulations and 
is calculated from the variance (AA; 2 ) in the Hamming distance distribution. Parameters are 
L = 8, a = 1, and P = 1. 




16 



the effects of specific replication on the error threshold. In our model, which 
accounts for the full quasispecies mutant distributions, replication rates depend 
on the Hamming distance between enzyme and substrate via an adjustable de- 
gree of specificity. Combining stochastic simulations, numerical solutions of 
reduced rate equations and analytical solutions to a binomial closure approx- 
imation, we could analyze the entire phase diagram and assess how mutation 
rate /i, selection strength (3 and specificity degree p influence the localization 
about a master sequence in order to preserve its information content. We found 
that for self-specific replication very weak specificity suffices for localization, 
whereas stronger specificity gives more tolerance against mutations but leads to 
bistability with the delocalized regime of random sequences. In particular, the 
binomial closure approximation permits to obtain analytical expressions for the 
bifurcation diagram and an upper bound /?* for the bistability regime, which can 
be evaluated for any specificity function f s (d). Apart from our special choice, 
Eq. (HJ , a mesa-shaped function would also be conceivable, in correspondence to 
fitness landscapes for transcription fa ctor binding allowing for some "fuzzincss" 
or neutrality in the binding sequence dGerland and HwaL 120021) . Preliminary re- 
sults indicate that in this case the binomial closure approximation gives at least 
qualitative agreement as well. While our approximation is not restricted to large 
L, this limit can probably be more systematically be describ ed using the max- 



imum principle employed pr eviously in quasispecies theory (jHermisson et al 



2002t ISaakian andHul . 120061 ). In the case of cross-specificity, we found that co- 



existence of subpopulations localized about complementary sequences is possible 
only if replication rates increase faster than linearly with Hamming distance. 

Although our model is based on idealizing assumptions, we can extrapolate 
from currently available experimental data to obtain rough quantitative esti- 
mates for the maximal length L c of the recognition or tag sequences that can 
be used by rcplicase enzymes to specifically and reliably discriminate appropri- 
ate and useless templates (and vice versa). Considering that non-enzymatic 
template-direc ted polymerization rates are on the order of several hours to 
days per base ( Acevedo and Orgel 1987 ; IWu and Orgel , 19921). while ribozyme- 



catalyzed polyme r izatio n gives rates in the hour range ([ Johnston et all 12001 



IZaher and Unraul 12007 ). we can estimate the ratio /3/a somewhere near 5-20 
if polymerization is the rate-limiting step. Assuming a self-specific enzymati c 



replicator with a mutation rate on the order of 3% as in ([Johnston et all 120011 ) 



we obtain a critical length L c = 11-15 for weak specificity p = 1, and a larger 
value L c = 18-33 for p — > oo, because stronger specificity constraints allow longer 
sequences due to better localization. These values are significantly smaller than 
the lengt hs of, e.g., t he 154- nucleotide specificity domain of Bacillus subtilis 
RNase P (jLillevl . 120051) or the tRNA-like s tructures supposed to act a s "genomic 
tags" for the replication of RNA viruses ( Weiner and Maizeld . Il987 ). Many of 
the nucleotides in these instances have a structural role, which makes them 



effectively redundant or neutral (|Kun et all 120051) , and only a minority is ac 



tually involved in recognition. Also, recent research indicates that "stalling" of 
polymerizatio n after mismatch incor poration might significantly reduce the er- 
ror threshold ( Raiamani et al. . 2010l ). Nevertheless, our result suggests that the 
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error threshold puts hard constraints on the information content of enzymatic 
replicators as well. 
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Appendix A. Self-specific replication 



Appendix A.l. Derivation of Eq. @ 

To obtain Eq. ©, wc compute the first moment J^k ^Xk of the reduced rate 



equations Eq. © under the assumption that Xk 
binomially distributed. This gives four terms 



x\ EE - a) 



L-k 



IS 



Ti 



jk 



T 2 = y^^km kj b jl x ] °x ] ° 



x b x b 



(A.l) 



The reduced mutation matrix rrikj and the catalytic matrix bij are given in 
Eqs. ([7]) and (jHJ). Note that although not immediately obvious, the catalytic 
matrix bij = bji is symmetric, because the binomial in the denominator of 
Eq. ([8]) normalizes it to the single sequence level. Wc further keep in mind that 
(?) = if k < or k > n if n and k are integer, so we do not need to keep track 
of the summation limits in the following calculations. 

While it is straightforward to find T\ = aL(a + /Lt(l — 2a)) and T3 = —aaL, 
we concentrate first on T4, which reads after performing the summation over k 



Ta 



L 



j + n 



f s (i-j+2n)a^{l-a) 



\2L-{i+j) 



(A.2) 



Replacing i> = i - j + 2n and rearranging (^) (J) (/j = (£) Q 

sum over j and n, and are left with 



Ta 



-PaL 



f s (i')(l - 2o(l - a)) L -' 1 ' (2o(l - a)) 1 



we can 



(A.3) 



-paLS{l - 2a(l - a)). 



The term T2 can after similar rearrangements be written as the product of two 
generalized Vandcrmonde matrices: 



Jfc 



1-A* 



(1 



a(l — a) 
a(l — a) 



./■' 



(A.4) 
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where the L Vandermonde matrix with parameters a, b, c, and d is defined as 



a b 
c d 



X] ( ,• p 



3 \ ~ L + e ~ i -jb j ~ e c i ~ e d e . 



(A.5) 



This allows us to use a nice multiplication identity (jRawlings and Szd . 120051 ) for 
these matrices: 

(A.6) 



a b 




e /" 




ae + bg 


af + bti 


c d 




g h 




ce + dg 


cf + dh 



which gives: 



r a =/>£(J)/.(0£* 



a 2 + fi(l - 2a) o(l-o)' 
(1 - a) 2 - p(l - 2a) o(l-o) 



L-k, 



=/3[L(a + fj,(l - 2o))5(l - 2o(l - a)) 

- (1 - 2a)a(l - a)(l - 2/z)5'(l - 2a(l - a))] . 
Adding up T x + T 2 + T 3 + T 4 = gives Eq. ©. 



(A.7) 

(A.8) 
(A.9) 



Appendix A. 2. Solutions of Eq. ^ 

This equation can be solved whenever S(x), defined in Eq. (|10|) . assumes a 
simple form. For f s (d) = (1 — d/L) p as in Eq. S{x) is a polynomial of order 
L, except for integer < p < L, where it is of order p. A few instances for 
L > 2 are given by 



S(x) 



l-(l-x) L , p = 0, 
x, P=l, 



L-l 2 

x", p = oo. 

To find solutions a(p) ^ 1/2, we write a; = a(l — a) and solve Eq. 

l/L 



(A.10) 



for 



x = < 



2 (2^) , 



2 

a+,3 



/3 l+2/i(L-l)' 



(A.ll) 



2L ' 



For the last case, we approximated (1 — 2x) J 



p = 
p = l 
p = oo. 

3 -2a;L and fx/(l — 2(1) fH p 

for the important asymptote L > 1 with pL fixed, and used Lambert's W- 
function. A more complicated expression is obtained for p — 2. The solution 
a = | (l ± VI ~ 4a; ) is then easily computed, and the error thresholds yU c follow 
from evaluating the condition x = 1/4 (for p = 0, 1,2), or from requiring the 
W-function to give a real result (for p = oo): 

/^ 2 -(L + l) ; p = Qj 



a+/3 



AaL+2[3(L+l) ' 
g 



P 



1, 



4aL+/3(L+3) ' 

J.' 

ea 



P = 2, 
p = oo. 



(A.12) 
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Finally, it is easy to evaluate the bifurcation diagram Eq. (fTTj) at a = 1/2 to 
get the critical mutation rate fl c = fi(l/2): 



/'c 



(3S'(l/2) 



20S'(l/2) + 4L(a + /3S(l/2))~ 
For our specificity function / s (d) = (1 — d/L) p , we obtain explicitly 



(A.13) 



f o-q+n P 

P 

4aL+2/3(L + l) ■ 

P 

AaL+fi(L+3) ' 
P 

I a2 L + 1 +ip ' 



p = 

p = l 

p = 2 
p = oo. 



(A.14) 



Note that fi c = fl c for p < 2. 

While most of these results can be evaluated also for a = 0, the case p — > oo 
is special. Here, Eq. ((9| gives a; = /z/(l— 2/i) if we again approximate (1 — 2x) L ps 
e~ 2:rL , hence the error threshold is fj, c = 1/6, independent of and L. 



Appendix B. Cross-specific replication 

Appendix B.l. Derivation of an equation for a 

To obtain an equation for the parameter a, we compute the second moment 
^2kk 2 ik of the reduced rate equations, Eq. ([6]), under the assumption that 
Xk = x% = ^ (fc) t -* (1 — a) L ~ k + a L ~ fe (l — a) k ] is a sum of two complementary 
binomials, because in this case the first moment vanishes by construction. The 
four terms T1-T4 are defined and evaluated analogously to Eq. (|A.1[) . and after 
some algebra we find: 



1 



Ti = -aL [L - 2(L - l)(a(l - a)(l - 2fx) 



T 2 = -0 

2 H 



T 3 



m(i-m))] 

L(L - 2{L - l)(a + /i(l - 2a))(l - a - /i(l - 2a)))C(l 
+ 2a(l - a)(l - 2a) 2 (l - 2^) 2 (L - 1)C"(1 - 2a(l - a)) 
+ 2(o(l - a)(l - 2a)(l - 2^)) 2 C"(1 - 2a(l - a)) 

- -<xL[L-2(L- l)o(l - a)] 

- -f3L[L - 2(L - l)a(l - a)] (7(1 - 2a(l - a)). 



(B.l) 
2a(l - a)) 

(B.2) 
(B.3) 
(B.4) 
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Here, we have defined C(x) = ±£ fc (£) [f c (k) + f c (L- k)]x L - k (l- x) k . Adding 
up Ti + T 2 + T 3 + T 4 = gives the condition 



(1 - 2a) 2 - n)L(L - 1) [a + /?C(1 - 2a(l - a))] 

- @a(l - a)(l - 2^) 2 [(L - 1)C*'(1 - 2a(l - a)) 

+ a(l -a)C"(l-2a(l- a))] I =0. (B.5) 



Most importantly, Eq. (|B.5[) reads for p = 1 in the specificity function / c (d) = 
(d/L)P: 

- (1 - 2afL{L - l)/i(l - + /3/2) = 0, (B.6) 
which has only the solution a = 1/2. 

Appendix B.2. Bifurcation diagram 

Solving Eq. (|B.5I) for /i gives the bifurcation diagram: 



1± 1 



4^a(l-a)[(L-l)C'(l-2a(l-a))+a(l-a)C"(l-2a(l-a))] 
L(L-l)[a+l3C{l-2a(l-a))\ 



-1/2 



(B.7) 

where we take the negative sign and the positive root to obtain values /x near 
zero (values near unity imply complementary replication and give equivalent 
results for cross-specific replication). 

We readily find that /x"(l/2) > if (3 > 0, which implies that the pitchfork 
bifurcation described through Eq. (|B.7I) is always subcritical. 



Appendix B.3. Solutions of Eq. (|B.5|) 

The auxiliary function C(x) can be evaluated for small integer p as in 
Eq. (| A. 101) . There is no solution to Eq. (|B.5|) except a = 1/2 for p = 1, and the 
expression for p = 2 is quite lengthy. For p ~ oo and a > 0, we get 



* = -(l-a) = M -^w(^e^ 



(Bi 



which is exactly the result for the self-specific case if we replace (3 — > (3/2. 
Accordingly, we get /j, c = W[/3/(2ea)]/(2L) for the error threshold. 

Observing that C"(l/2) = 0, the critical mutation rate jj c = ^(1/2) is given 

by 



(3C"{\/2) 



4L(L-l)(a + /3C(l/2)) 



-1/2' 



(B.9) 
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which reads explicitly 



p = l 

p = 2 (B.10) 
p = oo. 

We can simply take a = in most of the above expressions to investigate the 
case of zero non-enzymatic replication rate. In the limits p — > oo and L — > oo, 
we find that the bifurcation diagram /z(a) converges towards 

H(a) a(l - a), (B.ll) 

except for a region near a = 1/2, because /x(l/2) has to coincide with the exact 
value /Et c = (2 — \/2)/4. Because this region becomes infinitesimally small as 
L — > oo, we conclude that fi c — > 1/4 in this limit. 
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